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Abstract. A previous MHD theory for the density jump at the Earth’s bow shock, which 
assumed the Alfven ( Ma ) and sonic (M s ) Mach numbers are both >> 1, is reanalyzed 
and generalized. It is shown that the MHD jump equation can be analytically solved much 
more directly using perturbation theory, with the ordering determined by Ma and M s , and 
that the first-order perturbation solution is identical to the solution found in the earlier 
theory. The second-order perturbation solution is calculated, whereas the earlier approach 
cannot be used to obtain it. The second-order terms generally are important over most of 
the range of Ma and M s in the solar wind when the angle 0 between the normal to the bow 
shock and magnetic field is not close to 0° or 180° (the solutions are symmetric about 90°). 

This new perturbation solution is generally accurate under most solar wind conditions at 1 
AU, with the exception of low Mach numbers when 6 is close to 90°. In this exceptional 
case the new solution does not improve on the first-order solutions obtained earlier, and the 
predicted density ratio can vary by 10-20% from the exact numerical MHD solutions. For 
0 ~ 90° another perturbation solution is derived that predicts the density ratio much more 
accurately. This second solution is typically accurate for quasi-perpendicular conditions. 

Taken together, these two analytical solutions are generally accurate for the Earth’s bow 
shock, except in the rare circumstance that Ma < 2. MHD and gasdynamic simulations 
have produced empirical models in which the shock’s standoff distance a, is linearly related 
to the density jump ratio X at the subsolar point. Using an empirical relationship between 
a s and X obtained from MHD simulations, a s values predicted using the MHD solutions 
for X are compared with the predictions of phenomenological models commonly used for 
modeling observational data, and with the predictions of a modified phenomenological model 
proposed recently. The similarities and differences between these results are illustrated 
using plots of X and a s predicted for the Earth’s bow shock. The plots show that the 
new analytic solutions agree very well with the exact numerical MHD solutions and that 
these MHD solutions should replace the the corresponding phenomenological relations in 
comparisons with data. Furthermore, significant differences exist between the standoff 
distances predicted at low Ma using the MHD models versus those predicted by the new 
modified phenomenological model. These differences should be amenable to observational 
testing. 


1. Introduction 

Research on the position of the Earth’s bow shock, 
which has been pursued for the last 3 decades, has pre- 
dominately focused on gasdynamic and quasi-gasdynamic 
analysis [e.g., Spreiter et a/., 1966; Fairfield , 1971; 
Formisano et a/., 1971; Holzer and Slavin , 1978; Slavin 
et al , 1983; Farris et a/., 1991]. There are strong mo- 
tivations, however, for further developing an MHD the- 
ory of the shock’s location [Cairns and Grabbe } 1994], 
since the bow shock is a magnetosonic shock wave and 
the solar wind interacting with the Earth’s magneto- 
sphere is a magnetized plasma. Robust MHD theories 
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show dependence on the Alfven Mach number Ma , the 
sonic Mach number M $1 and the angle 9 between the 
magnetic field and shock normal vectors. In contrast, 
gasdynamic theory is a function of the lone Mach num- 
ber M g (phenomenological substitutions of Ma or M mt 
for M s are often made to model the Earth bow shock in 
variants on gasdynamic theory) and its validity is gen- 
erally limited to M\ » M, » 1. (Here M mt is the 
0-dependent magnetosonic Mach number.) Our focus is 
on developing an accurate model of the solar wind in- 
teraction for low Ma and M,, which requires an MHD 
analysis. 

An earlier effort at a theoretical analysis of the bow 
shock using MHD theory was pursued by Zhuang and 
Russell [1981]. They showed that the jump conditions 
lead to a complicated cubic equation for the density 
jump ratio as a function of the magnetic field compo- 
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nents and temperature in the solar wind. Their pa- 
per presents an involved method for solving that equa- 
tion for large Ma and M t in which the nonlinear alge- 
braic equation for the density jump ratio X = psw/pba 
is converted into 30 linear equations ( p sw and pf, 3 are 
the mass density of the upstream solar wind and im- 
mediately downstream of the bow shock, respectively). 
These equations were solved to obtain analytical solu- 
tions to first order in the Mach numbers, both of which 
were assumed to be >> 1. 

There is interest in the density jump ratio X for two 
reasons. One is intrinsic, since spacecraft experiments 
measure X directly at the bow shock. The other is its 
relationship to the standoff distance of the bow shock. 
Empirical analysis has shown the magnetosheath thick- 
ness (the bow shock standoff distance less the known 
distance to the magnetopause) is proportional to A, so 
knowing X allows predictions for the standoff distance 
[Spreiter et al } 1966; Cairns and Grabbe , 1994; Cairns 
and Lyon , 1995]. The difficulty with using the Zhuang- 
Russell analytical solution for comparison with data on 
the bow shock is that when Ma or M $ is sufficiently 
below 10, their solution generally underestimates X . 
Our analysis indicates that the underestimation can be 
significant for all 0 except for 9 relatively close to 0. 
This inaccuracy degrades the usefulness of the theory 
in modeling most shocks which are observed at these 
lower Mach numbers. 

Historically the trend in investigating the bow shock 
position has been to utilize gasdynamic results for com- 
parison with the observed shock standoff distances, al- 
though efforts were made to introduce some MHD con- 
cepts into the gasdynamic model. Spreiter and Rizzi 
[1974] present MHD solutions for solar wind magnetic 
fields aligned with the flow (0 = 0), in which case the 
MHD equations reduce to a set of equations directly 
analogous to the gasdynamic equations except for a 
modified equation of state. Spreiter and Stakara [1980] 
introduced a global quasi-gasdynamic model oriented 
toward its usability for comparison with data in which 
the flow is determined by solving the conservation equa- 
tions with the magnetic field terms neglected, and the 
magnetic field is then determined from the gasdynamic 
flow using the magnetic induction equation. 

Russell [1985] argued how the difficulty of comparing 
with MHD theory made the quasi-gasdynamic model 
more usable for data comparison and that the single 
Mach number in the quasi-gasdynamic theory (which 
is the sonic Mach number M s in the gasdynamic equa- 
tions) should be replaced with the magnetosonic Mach 
number Af m ,, as was done previously by Formisano 
et al. [1971] when modeling shock locations. This 
phenomenological model, in which M ms is substituted 
for M s in the gasdynamic equation obtained for a, by 
Spreiter et al. [1966], has often been used to analyze 
spacecraft data for the Earth bow shock [e.g. Farris et 
al. t 1991]. After comparing some observed bow shock 
crossings with gasdynamic theory, Slavin et al [1983] 


state that there appears to be no theoretical support 
for this substitution, and that the arguments advanced 
for the substitution appear valid only under conditions 
of strong IMF intensity or low Ma • Cairns and Grabbe 
[1994] demonstrate the lack of theoretical foundation 
for this substitution. Russell and Zhang [1992] and 
Cairns et al [1995] have shown that this phenomeno- 
logical model is observationally inadequate at low Ma 
and M ms < 3, thereby indicating a need for further 
developing the theory at these low Mach numbers. Re- 
cently, Farris and Russell [1994] proposed a modified 
phenomenological model for a s for low Ma • 

These developments and the continued reliance on 
phenomenological models for data comparison clearly 
indicate a need for further developing MHD theory, 
especially at low Mach numbers. Cairns and Grabbe 
[1994] considered the obvious MHD generalization of 
Spreiter et al's [1966] gasdynamic theory for a s : lin- 
early relating the MHD density jump ratio X (rather 
than the gasdynamic X) to the magnetosheath thick- 
ness. They then showed that the resulting theory for 
a s reduces to the phenomenological M mt model only 
under restricted circumstances which are generally not 
met in the solar wind. Their theory shows direct de- 
pendences on Ma , M S} 9 , and specific heat ratio *) 
that cannot be subsumed into M m , alone, and reduces 
to the gasdynamic theory when Ma » M t >> 1 
Subsequently, Cairns and Lyon [1995] presented three- 
dimensional ideal MHD simulations for 9 > 45°, which 
show that the magnetosheath thickness is indeed pro- 
portional to the MHD value of X for Ma > 1*5. This 
was utilized to extend the MHD model for a 9 developed 
by Cairns and Grabbe [1994]. 

The present paper markedly extends the Cairns ana 
Grabbe [1994] MHD analysis by obtaining analytical so- 
lutions of X for all 0, linking them with the MHD sim- 
ulations of Cairns and Lyon [1995], and comparing the 
predictions with the modified phenomenological model 
of Farris and Russell [1994]. In section 2 it will be 
shown that the MHD jump conditions can be analyt- 
ically solved, in a much more direct fashion than was 
done earlier, by using perturbation theory. The solu- 
tion given by Zhuang and Russell [1981] is the same as 
the first-order perturbative solution which is obtained 
much more directly. Furthermore, with perturbation 
theory the second-order (and even higher) terms in the 
solution are obtained, whereas they cannot be in the 
approach used by Zhuang and Russell. Two solutions 
accurate to second order are obtained for separate over- 
lapping ranges of 0: approximately 0° to 75° and 50° 
to 90° (with symmetry about 90°). The second- and 
higher-order terms are determined to be important for 
0 not close to 0° except when both Ma and M s are 
>> 1. Taken together, these two analytical solutions 
are generally accurate to within a few percent under 
most conditions that exist in the solar wind at 1 AU. 
Section 3 compares the low-Mach number predictions 
from both the new analytic theory, and from magne- 
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tosonic phenomenological models that have been used 
in spacecraft data analysis, against the precise MHD re- 
sults. These comparisons are made between the shock 
standoff distances predicted by selected empirical mod- 
els. Section 4 summarizes the major conclusions of this 
paper. 

2 • Density and Velocity Jump 

The Rankine-Hugoniot relations for the conservation 
of mass, energy, momentum, and charge density, along 
with Maxwell’s equations, constrain the changes in fluid 
and electromagnetic variables across the shock. Zhuang 
and Russell [1981] used these relations to derive a cu- 
bic equation for the upstream to downstream density 
ratio X = p 9W /Pbi = Ub s /U 9W (sw denotes the solar 
wind upstream and bs the downstream side of the bow 
shock, while the symbols U denote fluid speeds). This 
equation has some differences from a similar cubic de- 
rived for shock waves using the Vlasov (rather than a 
fluid) approach of Tidman and Krall [1971]. The Vlasov 
derivation assumed Maxwellian electron and ion distri- 
butions, which are probably not valid near the Earth’s 
bow shock. The fluid equation can be expressed in 
the following streamlined form [see Cairns and Grabbe , 
1994]: 

X 3 + C a X 2 + C b X + C c = 0 (1) 


cally solve (1) for the bow shock under the restriction 
that Ma and M, » 1. Unfortunately, that approach 
ties the solution up into a Gordian knot. To unravel 
this knot, they broke the cubic solution into 30 linear 
equations, which were solved simultaneously for large 
Ma and M t . 

A better approach for finding analytical solutions 
avoids using the cubic formula altogether. Since the two 
Mach numbers are both > 1 for Earth’s bow shock, a 
direct approach is to solve the equation using perturba- 
tion theory [e.g., Nay f eh , 1993], with the perturbation 
parameters being 1/Af^ and 1 /M 2 . These are two dif- 
ferent perturbation parameters, but they will be taken 
as the same order e in our analysis since they are ap- 
proximately of the same order in the solar wind. We 
will denote the order given by these two parameters as 
e. Thus the ordered solution we are seeking takes the 
form: 

X = x 0 + exi + £ 2 X 2 + ... (5) 

The coefficient terms in (1) are similarly ordered into 
powers of e: 

C a — - R a — tS a (6) 

C b = eS b + e 2 T b (7) 

C c = —e 2 T c — e 3 U c (8) 

where the new variables are 


where the coefficients are 
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where Ma = U 9W /Ua and M s = U sw /c t are the Alfven 
and sonic Mach numbers, 9 is the angle between the 
shock normal and the solar wind magnetic field B , and 
7 is the adiabatic constant. Here U sw is the solar wind 
speed, Ua is the Alfven speed, and c s the ion acoustic 
speed. The analysis by Cairns and Grabbe [1994] will be 
extended by finding analytical solutions over the entire 
range of 9 . 

While (1) has three solutions for any set of plasma 
parameters, only one corresponds to a solution for the 
Earth’s bow shock. For 9 close to 0°, the other two so- 
lutions correspond to “switch-on” shocks ( X ~ 1/M a), 
while for 9 in the vicinity of 90° the other two solutions 
represent an unphysical “negative- density” shock (mass 
density of opposite sign on the two sides of the shock) 
and an unphysical “reverse” shock (X < 0.25). Zhuang 
and Russell [1981] used the cubic formula to analyti- 
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From these orderings the perturbation equations for 
0(e n ) can be written in a straightforward manner. The 
0(e°) equation just gives an expression for xq, which 
can then be substituted into the 0(e) equation to alge- 
braically solve for x\\ 


x 0 = R a (15) 

x\ — S a - Sb/R a (lb) 

Combining Xo and Xi then gives a solution for X to first 
order in e: 

y — (7 ~ 1 2 r 1 1 sin2g 1 / 17 i 

(7+1) (7+l)] l M2 + ( T -l)M2 J 1 } 

This solution is precisely the same as that obtained by 
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Zhuang and Russell [1981, equation 43] through their 
analysis with 30 linear equations. 

In addition to the more direct route in obtaining the 
first-order solution, the advantage of the perturbation 
technique is that the solution can be extended to the 
next higher order in e without difficulty. The second- 
order equation is 

3a? 0 (^i + x 0 x 2 ) - [R a (x\ 4- 2x 0 x 2 ) (18) 

-f2xo£iS<i] + (xo Tft + xi Sb) — T c — 0 

Solving this equation, the solution for X to second order 
becomes 

X _ (mi) ■ 2 r 1 ■ sin2g I / 19 x 

" (7 + l) + (7+l) l A /. 2 + (7-l)M3 J ^ 

2 sin 2 0 , 7 cos 2 0 . 1 

+ (7-l) 2 Ml ! (7-l jM\ +{y ~ 

1 n 

+ (7-l)Mj J) 

Note that the second-order term in (19) disappears if 
9 = 90° and 7 is taken as 2. This characteristic behav- 
ior may have played a role in the conclusion that 7 = 2 
which Zhuang and Russell obtained by comparing data 
with their first-order theoretical results. Their data set 
could well have contained a preponderance of quasi- 
perpendicular shocks. Thus the limitations of their 
first-order theory may have created an inherent bias in 
the data comparison toward an extrapolated value for 7 
which minimizes the second-order corrective term, that 
is, toward 7 = 2 for cases with 9 ~ 90° . This was very 
close to the average value actually determined in their 
analysis. Later data analyses [e.g., Farris et a/., 1991] 
have concluded that the value of 7 is actually closer to 
5/3. 

In fact, when 9 is close to 90° the foregoing per- 
turbation expansion becomes relatively inaccurate be- 
cause the contribution of the second-order terms be- 
comes small (the 7 — 2 terms become predominant), and 
thus third- and higher-order terms in the perturbation 
expansion become important except when both Mach 
numbers are >> 1. In this case, however, the solution 
can be replaced with a better one. In (1) the coefficient 
C c and two terms in C& become of third or higher order 
in a modified perturbation ordering because the cos 2 6 
factor becomes very small (going to 0 at 9 = 90°). By 
neglecting these third- or- higher order terms the equa- 
tion for X becomes quadratic and is easily solved: 

x = ^{ 1 + ( 20 ) 

\fi - 4[(7 + 7 cos 2 0 - 2)/(7 + 1 

where C a is just the coefficient given by ( 2 ). This so- 
lution is accurate to second or higher order for 9 away 
from 90° (as far down as 50° — 60° ), and collapses to 
the solution by Cairns and Grabbe [1994] as 9 — > 90°. 


Taken together, (19) and (20) constitute reasonably 
simple analytical relationships derived from MHD the- 
ory for the density jump across the bow shock. These 
expressions can be easily compared with data. Note 
that these two solutions have a clear dependence on 
four variables: Ma } Af, , 0, and 7 , as expected theoret- 
ically [Cairns and Grabbe , 1994]. This variable depen- 
dence is missing in previous theoretical models for the 
bow shock location that use expressions for the den- 
sity jump based on gas dynamics. Such models have 
previously predominated comparisons with data on the 
Earth’s bow shock location; the dependence on 4 vari- 
ables is only present in the Zhuang- Russell MHD model 
and in our recent MHD work [ Cairns and Grabbe , 1994; 
Cairns and Lyon , 1995; and this paper]. 

Predictions for X as a function of Ma are compared 
in Figures 1-5 for the analytical solution from either 
(19) or (20), the exact numerical MHD solution, and 
the first-order solution of Zhuang and Russell [1981]. 
It should be noted that Zhuang and Russell derived 
their results under the assumptions that Ma and M, 
are >> 1. Therefore curves of their results for the low 
Mach numbers are used here only for the purpose of 
comparison with the predictions of (19) and (20), whose 
validity definitely does extend to lower Mach numbers. 

In Figures 1-3 the analytic solutions for X as a func- 
tion of Ma are compared for four angles between the 
magnetic field direction and the shock normal ( 9 = 
0°, 30°, 60°, 90°) and for three different sonic Mach num- 
bers (M, = 8,4,2). Figures 1 and 2 span the range of 
M a generally observed in the solar wind at 1 AU [e.g., 
Fairfield , 1971, Figure 3]. Figure 3 is for a value of M $ 
just below the most extreme cases observed in the solar 



Figure 1. Comparison of the density ratio X — 
p$w / Pb$ predicted by the second-order perturbation 
MHD solutions (19) and (20) (dashed lines) with the 
exact numerical MHD solution (solid lines) and the 
Zhuang- Russell solution (dotted lines), as a function of 
Alfven Mach number Ma and selected 9 (angle between 
the bow shock normal and the IMF) for sonic Mach 
number M 8 = 8 . The adiabatic constant 7 = 5/3 is 
used in all plots. 
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wind and is used to test how far the theory extends. 
All three figures show that the second-order solutions 
(19) and (20) predict the exact numerical results more 
accurately than the first-order Zhuang-Russell solution 
for all angles except for 0 = 0°, where both the new an- 
alytical solution and the Zhuang-Russell solution agree 
precisely with the exact numerical solution. 

Figures 1-3 show that the Zhuang-Russell result for 
the low Mach-numbers becomes increasingly inaccurate 
as M\ decreases for a given 9 , and that the their solu- 
tion underestimates the exact solution for finite 9. The 
Zhuang-Russell solution shows a noticeable error for all 
three nonzero angles in each of Figures 1-3. The ana- 
lytic solutions contained in (19) and (20) show relatively 
small errors at intermediate 9 (30° and 60°) when Ma 
goes below 2-3. 

Figures 4 and 5 show how these predictions for X 
vary with 0. In Figure 4 the value of X given by (19) is 
plotted as a function of 9 (ranging from 0° to 70°) for 
M s = 4 and selected values of Ma • These plots show 
that the error in the analytic theory becomes significant 
at most angles when Ma — 2, although they are gener- 

o.e 
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Figure 3. Same as Figure 1 for Af, = 2. 




e 

Figure 4. Comparison of the density ratio X predicted 
by (19) with the exact numerical MHD solution and 
the Zhang- Russell solution, as a function of 9 ranging 
0° — 70° for M s = 4 and selected small values of Ma • 

ally negligible for higher In all cases (19) is more 
accurate the Zhuang-Russell result. 

In Figure 5 the value of X given by (20) is plotted as 
a function of 6 (ranging from 60° to 90°) for M 3 = 4 
and selected values of Figure 5 shows that (20) 
is very accurate (within about 1-2%) for Ma >1.5 for 
the approximate range of 9 — 75° — 90° and only shows 
more significant error in the range 9 = 60° —75° for the 
case that Ma = 1.5. In all cases the first-order result 
of Zhang and Russell is less accurate. 

3. Analytical Solutions and 
Phenomenological Models 

Equations (19) and (20), which were derived analyti- 
cally from MHD theory, constitute useful analytical so- 
lutions that both complement the analysis by Cairns 
and Grabbe [1994] and can be easily compared with 
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Figure 5. Comparison of the density ratio X predicted 
by (20) with the exact numerical MHD solution and 
the Zhang-Russell solution, as a function of 9 ranging 
60° — 90° for M a = 4 and selected small values of Ma . 
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shock crossing data. As discussed in section 1 gasdy- 
namic [Stiff, 1962; Spreiter et al., 1966] and MHD simu- 
lations for 0 > 45° [ Cairns and Lyon , 1995] show empir- 
ically that the standoff distance a, is linearly dependent 
on the density jump X. (Note that a B is defined here 
as the distance between Earth’s center and the shock’s 
nose; in the aerodynamic literature ’’standoff distance” 
is defined as the distance of the obstacle to the shock, 
so an alternative would be to call a $ the “geocentric 
shock distance.”) This section examines and compares 
the predictions of models for a,, that are constructed 
from MHD and gasdynamic solutions for X combined 
with the empirical relationships between a, and X, for 
parameter ranges that can occur at the Earth’s bow 
shock. 

These empirical relationships between a 9 and X = 
psw/pbs take the following form [Seiff, 1962; Spreiter 
et al. , 1966; Farris and Russell , 1994; Cairns and Lyon , 
1995]: 


— = j + kX (21) 

a mp 

Here a mp is the distance from Earth to the nose of the 
magnetopause. It should be emphasized that although 
the analytical solutions found in section 2 for the density 
ratio X apply locally everywhere on the bow shock, 
these empirical models are restricted to the subsolar 
region (the nose) of the shock. In the more general 
case, a, depends on (at least) five parameters: the four 
parameters (Ma> M s , 0 , and 7) that X depends on, as 
well as the magnetopause shape [Fpneifer and Stahara, 
1980; Slavin et al., 1983]. 

For the gasdynamic empirical relation found by Seiff 
[1962] and further developed by Spreiter et al. [1966] 
j = 1 and Jfc = 1.1, where the value of k depends on the 
obstacle shape. This model has been the prevailing one 
used for published comparisons with spacecraft data for 
many years. Solving the gasdynamic equations for the 
density jump X yields 


(7 - 1 )M[ ± 2 

(7 + 1)3/2 


( 22 ) 


The phenomenological model of Formisano et al. [1971] 
and Farris et al. [1991] substituted M mt for M, . See 
Cairns and Grabbe [1994] for arguments against this 
procedure. 

Two alternative empirical models have been devel- 
oped recently, principally due to difficulties account- 
ing for distant shock locations at low M a and M ms 
[Russell and Zhang , 1992; Cairns et al., 1995]. In the 
model presented by Farris and Russell [1994] the value 
for k is modified from the Seiff- Spreiter form at lower 
Mach numbers by it = 1.1 — 1) while j stays 

1. In the model developed from MHD simulations by 
Cairns and Lyon [1995], j — 0.4 and k = 3.4 for quasi- 
perpendicular flows with M s ~ 8 and Af^ > 1.5. This 
model explains their simulation results with excellent 
accuracy. Both the alternative models were developed 


with the goal of extending the relationship (21) to lower 
Mach numbers and exhibit a greater dependence on the 
density ratio at these low Mach numbers. It should 
be noted that for larger Mach numbers (> 5 — 10), 
psw/pit ~ 0.25 and all three empirical models are in 
approximate agreement. 

A note regarding the importance of 0 for these em- 
pirical models is in order. Cairns and Lyon [1995] point 
out that a t depends strongly on 0, as evidenced by the 
qualitatively opposite variations of a, with decreasing 
Ma at fixed M, for their MHD simulations (valid for 
0 > 45°) compared to the field- aligned (0 = 0°) MHD 
simulations of Spreiter and Rizzi [1974]; that is, a, is 
predicted to increase with decreasing Ma for 9 > 45° 
but to decrease with decreasing Ma for 0 = 0°. The 
phenomenological models also predict variations in a t 
opposite to the Spreiter-Rizzi field-aligned MHD simu- 
lations, and this characteristic also argues against the 
validity of the phenomenological models near 0 = 0°. 
It is presently unknown how large 0 must be for mod- 
els based on (21) to hold, but the empirical models are 
likely relevant for 0 > 20°. The recent observational 
results of Peredo et al. [1995], who find statistically 
that a, decreases with decreasing Ma (> 2) but do not 
specify 0, make the dependence of a, on 0 particularly 
topical. 

In Figure 6 comparisons are made between the pre- 
dictions for a,/a mp , based on the Cairns-Lyon empir- 
ical values of j and Jt, for various 0 (0°, 30°, 60°, 90°) 
at M t = 8. The solid line shows the exact numerical 



Figure 6. Comparison between the second-order an- 
alytical MHD solution (dashed lines), the numerical 
MHD solution (solid lines) and the gasdynamic solu- 
tion with a phenomenological substitution of the mag- 
netosonic Mach number Af mj for the sonic Mach num- 
ber M s (dotted line). Here 7 = 5/3 and M t = 8. 
The axis on the left-hand side is a 9 /a mp predicted for 
both models using the empirical relation between X and 
found in Cairns and Lyon’s [1994] simulations. 
It is expected that this empirical relation is wrong for 
the 0 = 0° case (from the Spreiter and Rizzi [1974] stud- 
ies at 0 = 0°) so the value of X from which the standoff 
distance is determined is shown on the right-hand side. 
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m a 

Figure 7. Same as Fig. 6 for the case that M a = 4. 


MHD solutions, while the dashed line uses the analyt- 
ical solutions (19) and (20) for X. The dotted line 
shows the prediction for the phenomenological M ms 
model of Formisano et ai [1971] and Farris and Rus- 
sell [1991] (substitute M ms for M s in (22) and substi- 
tute the result into (21)) with the modification j = 0.4 
and k = 3.4. Note that Cairns and Lyon [1995] show 
that, on the other hand, retaining the SeifF-Spreiter re- 
lationship j = 1 and Jr = 1.1 yields predictions for the 
phenomenological model that are widely different from 
their M a < 5 simulation results. Since the empirical 
models are not valid for small Q, that is, £ 20° (see pre- 
vious paragraph), a scale corresponding to X is given 
on the right-hand side. 

Similar comparisons are made in Figure 7 for M, = 4. 
As expected from Figures 1-3 the predictions from (19) 
and (20) differ from the exact MHD results only for 
Ma ~ 2 — 3 at intermediate 9 (30° and 60°) but agree 
very well down to M A = 1.5 at 9 = 0° and 90°. Predic- 
tions based on the phenomenological model, however, 
show much larger deviations from the exact MHD solu- 
tion at 9 = 0° and 9 = 90° . For the intermediate angles 
the phenomenological M ms model exhibits more mod- 
est deviations from the numerical MHD solutions but 
still generally show larger differences than the analytic 
solutions do. As 9 increases toward 90°, the predictions 
of the phenomenological model are initially larger but 
eventually smaller than the MHD results for all Ma- 
These differences appear suitable for experimental test- 
ing. 

For the case 9 — 0° , where both the analytic solution 
and the earlier Zhuang-Russell solution agree precisely 
with the exact MHD solution, Figures 6 and 7 show that 
the magnetosonic phenomenological model (replace M a 
in (22) by M mt ) exhibits a greater than 25% over- 
shoot of the MHD solution for X at Ma = 2. Thus 
the phenomenological model exhibits a dramatic fail- 
ure at low 9 , as was pointed out by Cairns and Grabbe 
[1994, Figure 5]. Note ironically that the size of this 
error increases as M 3 increases. The phenomenological 


”Alfven” model (substitute Ma in (22) for Af,), which 
was sometimes compared with data in the 1970s [e.g., 
Fairfield i 1971; Formisano et a/., 1971], corresponds to 
the magnetosonic form at 9 = 0° extended to model all 
9. Note again that the predicted a $ scale is not reliable 
for any of the models at 0 = 0° since the empirical rela- 
tion (21) breaks down there. While combining the MHD 
expression for X with the empirical relation (21) leads 
to the MHD models predicting standoff distances that 
are independent of Ma, the work by Spreiter and Rizzi 
[1974] for 9 = 0° predicts that a, actually decreases 
with decreasing Ma ; the behavior found by Spreiter 
and Rizzi is exactly opposite to that predicted for the 
magnetosonic and Alfven phenomenological models. 

For the case 9 = 30°, the phenomenological model 
exhibits a smaller error than at 9 = 0°. It reaches a 
maximum overshoot of almost 10% at around Ma = 3 
for the M g — 8 case in Figure 6, but this maximum 
falls dramatically as M s decreases in Figure 7. The 
analytic solution is much more accurate, and only shows 
> 1% deviation as Ma becomes close to 2. For the case 
9 = 60°, the phenomenological model exhibits a few 
percent undershoot. The analytic model is very close to 
the exact numerical solution for Ma > 3 but exhibits a 
smaller overshoot at Ma < 3. 

For the 9 — 90° case the analytic solution is equal to 
the exact solution, but the magnetosonic model exhibits 
a sizeable error. While the error in the phenomenologi- 
cal model is only a few percent in Figures 6 and 7, there 
is no error in the analytic model for 9 = 90° and almost 
none for 9 close to 90° . Thus even though one can argue 
that the magnetosonic model is a reasonably good ap- 
proximation near 9 = 90°, the analytic MHD solution 
is clearly preferable because it is a much more accurate 
model and because it is based on MHD plasma theory 
rather than phenomenology. 



Figure 8. Comparison between the second-order an- 
alytical MHD solution (dashed lines), the numerical 
MHD solution (solid lines) and new phenomenological 
model of Farris and Russell [1994] (dotted line) using 
different empirical models. Here 7 = 5/3 and M s = 8. 
The axis on the left-hand side is the predicted ratio 

a s / a mp • 
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Figure 9. Same as Figure 8 for M s — 4. 


In Figs. 8 and 9 comparisons of the MHD analyti- 
cal solutions are made with the new phenomenological 
model of Farris and Russell [1994] for 9 = 0° , 30°, 60° , 
and 90°, for M A = 8 and M A = 4, respectively. Once 
again the analytical and numerical MHD solutions for 
X are combined with Cairns and Lyon' s [1995] empir- 
ical relation to predict a$/a mp . No modifications are 
made to the Farris and Russell model since it is in- 
tended to account for shock locations at low M A and 
M mt without change. It is obvious that the two models 
make very different predictions for all 9 and M A <10. 
These differences are clearly amenable for observational 
testing. It should be noted that in general the rela- 
tive differences between the MHD model and the phe- 
nomenological model are sizeably larger for M A £ 1.5 
at 0 = 30°, 60° and 90° than they were in Figures 
6 and 7 (which have identical j and k in (22)). The 
predictions of the Farris and Russell model generally 
undershoots those of the MHD model. The increased 
difference between Figures 6 and 7 and Figures 8 and 9 
indicate a major discrepancy between the predictions of 
the MHD model and the phenomenological model. Fig- 
ure 9 indicates that the discrepancy increases in size as 
M t gets smaller, and that the two models now almost 
never agree for all 9 . 


except when 9 is very close to 90° while M g or M A is 
relatively small. For 9 ~ 90° the solution (20) much 
more accurately predicts the value for the density ratio. 
This second solution is quasi-perpendicular in nature 
and is generally also of second order in perturbation, of 
higher order as 9 approaches 90°, and becomes exact at 
9 = 90°. 

Taken together, (19) and (20) constitute reasonably 
simple analytical relationships derived from MHD the- 
ory for the density jump across the bow shock. When 
combined with appropriate models at the subsolar re- 
gion (bow shock nose) for the relation between X and 
the bow shocks standoff ratio a 9 /a mp , they can be eas- 
ily used for comparison with data. Two recent em- 
pirical models were used for comparisons and definite 
differences were found in their predictions. Equations 
(19) and ( 20 ) more accurately describe exact theoreti- 
cal MHD results than the phenomenological models for 
X that have been used (in conjunction with the Seiff- 
Spreiter empirical model) for comparison with space- 
craft data for many years. Those two equations, when 
combined with the most appropriate empirical model, 
constitute decidedly better replacements to use for that 
purpose. 

A goal for future research in this area is to use MHD 
theory to develop a fully analytic model for the re- 
lationship between X and a t /a mp to replace the em- 
pirical models. Only when that is accomplished will 
a completely theoretical MHD model be available to 
predict the bow shock’s location. Developing a full 
MHD description for the solar wind flow at the Earth’s 
bow shock, which will require a self-consistent three- 
dimensional determination of the shock location, shape, 
and jump conditions as a function of M A , M 9 , 7 , and 
9 , should then be possible. 
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